
#### Preliminaries ####

R.version$version.string

# clear workspace
rm(list=ls())

# load packages
library(here)
library(tidyverse)
library(rworldmap) # maps
library(magrittr)
library(reshape2)
library(ggpubr)
library(leaflet)
library(plotly)
library(gapminder)
library(countrycode)
library(ggrepel)
library(ggthemes)
library(modelsummary)
library(fixest)
library(marginaleffects)
library(estimatr)
library(brglm)
library(GJRM)
library(lmtest)
library(sandwich)
library(zoo)



ggplot_theme <- theme(axis.text=element_text(size=16),
                      axis.title=element_text(size=16,face="bold"), 
                      title =element_text(size=16, face='bold'))

legend_size_theme <-   theme(legend.key.size = unit(1, 'cm'), #change legend key size
                             legend.key.height = unit(1, 'cm'), #change legend key height
                             legend.key.width = unit(1, 'cm'), #change legend key width
                             legend.title = element_text(size=16), #change legend title font size
                             legend.text = element_text(size=16)) #change legend text font size 


#### set the working directory, if not using the Rproject #####

#---------------------------------------------------#
#--------------------load data ---------------------#
#---------------------------------------------------#

# 1. Load V-Dem data (version 14)

vdem_full_v14 <- readRDS("data/vdem14/V-Dem-CY-Full+Others-v14.rds")

vdem_subset <- vdem_full_v14 %>%
  dplyr::select(country_name, year, country_text_id, COWcode, country_id, starts_with("v2x_polyarchy"),
                starts_with("v2xca_academ"), starts_with("v2cafres"), starts_with("v2cafexch"), 
                starts_with("v2cainsaut"), starts_with("v2casurv"), 
                starts_with("v2clacfree"), e_pop, e_pop_sd, e_gdppc, e_gdppc_sd,  starts_with("v2x_liberal"), 
                v2x_regime, e_regionpol_6C, e_regionpol, v2canuni, v2caprotac, v2x_veracc, v2x_horacc, v2x_diagacc, 
                v2lgbicam, v2juhcind, v2jureview, v2exl_legitideol, v2exl_legitperf, v2exl_legitratio, 
                v2cacamps,e_boix_regime, v2petersch ) %>%
  filter(year >=1899)

episode_data <- read.csv("data/episode_data/episode_with_uncertainty_interval.csv")

episode_data <- episode_data %>%
  dplyr::select(-c(X, starts_with("v2xca_academ"), country_id, country_name))

summary(episode_data)
length(table(episode_data$decline_episode_id)) # 198 episodes

ert_data <- read.csv("data/ert_v14.csv")

ert_data <- ert_data %>%
  dplyr::select(-c(X, starts_with("v2x_polyarchy"), v2x_regime, country_id, country_name))

bmr_data <- read.csv("data/bmr/bmr_data_extended.csv")

#---------------------------------------------------#
#--------------------merge data --------------------#
#---------------------------------------------------#

vdem_subset <- vdem_subset %>%
  left_join(episode_data, by = c("country_text_id", "year")) %>%
  left_join(ert_data, by = c("country_text_id", "year")) %>%
  left_join(bmr_data, by = c("COWcode" ="ccode", "year")) %>%
  group_by(country_text_id) %>%
  mutate(v2xca_academ_lag = lag(v2xca_academ)) 

vdem_subset <- vdem_subset %>%
  distinct(country_text_id, year, .keep_all = T)

#---------------------------------------------------#
#----Autocratization Episodes Lags Construction ----#
#---------------------------------------------------#


vdem_subset <- vdem_subset %>%
  group_by(country_text_id) %>%
  mutate(aut_ep_lag_1 = lag(aut_ep), 
         aut_ep_lag_2 = lag(aut_ep, 2),
         aut_ep_lag_3 = lag(aut_ep, 3),
         aut_ep_lag_4 = lag(aut_ep, 4),
         aut_ep_lag_5 = lag(aut_ep, 5),
         aut_ep_lag_1_3 = ifelse(aut_ep_lag_1 == 1 | aut_ep_lag_2 == 1 | aut_ep_lag_3 == 1, 1, 0), 
         aut_ep_lag_1_5 = ifelse(aut_ep_lag_1 == 1 | aut_ep_lag_2 == 1 | aut_ep_lag_3 == 1 | aut_ep_lag_4 == 1 | aut_ep_lag_5 == 1, 1, 0)) %>%
  mutate(aut_ep_lag_1_3 = ifelse(aut_ep_lag_1_3 == 0 & aut_ep == 1, 1, aut_ep_lag_1_3), 
         aut_ep_lag_1_5 = ifelse(aut_ep_lag_1_5 == 0 & aut_ep == 1, 1, aut_ep_lag_1_5)) 

summary(vdem_subset$e_boix_regime)
summary(vdem_subset$democracy)

test_data <- vdem_subset %>%
  filter(year >=2019) %>%
  dplyr::select(country_name, year, e_boix_regime, democracy) %>%
  filter(is.na(democracy))

vdem_subset <- vdem_subset %>%
  mutate(democracy == ifelse(country_name == "Sudan" & year %in% c(2019,2020), 0, democracy), 
         democracy == ifelse(country_name == "Sudan" & year %in% c(2019,2020), 0, democracy))

summary(test_data)


#---------------------------------------------------#
#------------------Table A1 SA  --------------------#
#---------------------------------------------------#

#### Table A1 ####

democracies_afi_decline <- vdem_subset %>%
  filter(year >=2000 & !is.na(decline_episode_id)) 

democracies_afi_decline <- democracies_afi_decline %>%
  group_by(decline_episode_id) %>%
  filter(any(v2x_regime >=2))

democracies_afi_decline <- democracies_afi_decline %>%
  dplyr::select(country_text_id, year, decline_episode_id, v2x_regime, aut_ep_id, everything())

country_list <- unique(democracies_afi_decline$country_name) # 20 cases of AFI decline in demcoracies

democracies_afi_decline <- vdem_subset %>%
  filter(year >= 1998) %>%
  filter(country_name %in% country_list) %>%
  group_by(decline_episode_id) %>%
  summarize(country_name = first(country_name), 
            start_year = first(year), 
            end_year = last(year), 
            AFI_decline = first(v2xca_academ_lag) - last(v2xca_academ), 
            v2x_regime_start = first(v2x_regime), 
            v2x_regime_last = last(v2x_regime), 
            aut_ep = any(aut_ep ==1)) %>%
  drop_na(decline_episode_id) %>%
  dplyr::select(-decline_episode_id) %>%
  mutate(v2x_regime_start = case_when(v2x_regime_start == 0 ~ "Closed Autocracy", 
                                      v2x_regime_start == 1 ~ "Electoral Autocracy", 
                                      v2x_regime_start == 2 ~ "Electoral Democracy",
                                      v2x_regime_start == 3 ~ "Liberal Democracy"), 
         v2x_regime_last = case_when(v2x_regime_last == 0 ~ "Closed Autocracy", 
                                     v2x_regime_last == 1 ~ "Electoral Autocracy", 
                                     v2x_regime_last == 2 ~ "Electoral Democracy",
                                     v2x_regime_last == 3 ~ "Liberal Democracy"))

democracies_afi_decline$start_year <- as.character(democracies_afi_decline$start_year)
democracies_afi_decline$end_year <- as.character(democracies_afi_decline$end_year)


democracies_afi_decline <- democracies_afi_decline %>%
  arrange(aut_ep, country_name)

datasummary_df(democracies_afi_decline, output = "outputs/Table_A1.docx")

#---------------------------------------------------#
#------------------Table A2 SA  --------------------#
#---------------------------------------------------#

## not constructed by R-code but manually by literature review 


################################################################################
################################################################################

#---------------------------------------------------#
#------------------ Figure 1 -----------------------#
#---------------------------------------------------#

figure_1_data <- vdem_subset %>%
  filter(country_name %in% country_list) %>%
  filter(year >=1995) %>%
  mutate(linecol_aut = case_when(aut_ep == 1 ~ 1,
                                 T ~ 0),
         linecol_dem = case_when(dem_ep == 1 ~ 1,
                                 T ~ 0),
         linecol_decline = case_when(decline_episode == 1  ~ 1,
                                     T ~ 0),
         linecol_growth = case_when(increase_episode == 1 ~ 1,
                                    T ~ 0))

# define breaks and labels for x-axis for all plots
breaks  <-  seq(1995, 2025, 5)
labels <- ifelse(breaks %in% seq(1995, 2025, 5), breaks, "")
tick.sizes <- ifelse(breaks %in% seq(1995, 2020, 10), .75, .5)

ggplot(data = figure_1_data, aes(x = year )) +
  geom_line(aes(x = year, y = v2xca_academ, color = factor(linecol_decline), group=NA,
                linetype = "AFI"), alpha = .70, size = 1.1) +
  geom_line( aes(x = year, y = v2x_polyarchy,group=NA,
                linetype = "Polyarchy"), alpha = .70, size = 1.1, color = "#00AFBB") +
  geom_ribbon(aes(x = year, ymin = v2xca_academ_codelow,
                  ymax = v2xca_academ_codehigh), fill = "grey70", alpha = .45) +
  geom_ribbon(aes(x = year, ymin =v2x_polyarchy_codelow,
                  ymax = v2x_polyarchy_codehigh), fill = "#00AFBB", alpha = .45) +
  coord_cartesian(xlim = c(1995,2025), ylim = c(0,1), clip = "off") +
  scale_y_continuous("Academic Freedom Index", limits = c(0,1),
                     sec.axis = sec_axis(~ ., name = "Electoral Democracy Score"),
                     expand = expansion(mult = c(0, 0))) +
  xlab("") + theme_clean() + 
  facet_wrap(~country_name, nrow = 4, ncol = 5) +
  scale_x_continuous(breaks = breaks,
                     labels = labels) +
  scale_linetype_manual(name = "V-Dem Measure", values = c("solid", "dotted")) +
  scale_color_manual(name = "Episode", breaks = c(0,1),
                     labels = c("AFI","AFI decline Episode"),  values = c("black","#d7191c")) +
  theme(legend.position = "bottom")

ggsave("outputs/Figure_1.png", dpi = 600, width = 14, height = 10)

################################################################################
################################################################################

#---------------------------------------------------#
#------------------ Table 1 ------------------------#
#---------------------------------------------------#

democracies_afi_decline <- democracies_afi_decline %>%
  mutate(short_episode = ifelse((as.numeric(end_year) - as.numeric(start_year)) >= 4, 1, 0)) %>%
  group_by(short_episode) %>%
  summarise(n = n())


################################################################################
################################################################################

#---------------------------------------------------#
#------------------ Table 2 ------------------------#
#---------------------------------------------------#

democracies_autocratiaztion_cases <- vdem_subset %>%
  filter(year >=2000 & !is.na(aut_ep_id)) 

democracies_autocratiaztion_cases <- democracies_autocratiaztion_cases %>%
  group_by(aut_ep_id) %>%
  filter(any(v2x_regime >=2))

democracies_autocratiaztion_cases <- democracies_autocratiaztion_cases %>%
  dplyr::select(country_text_id, year, aut_ep_id, v2x_regime, decline_episode_id, everything())

country_list <- unique(democracies_autocratiaztion_cases$country_name) # 51 cases of autocratization 

democracies_autocratiaztion_cases <- vdem_subset %>%
  filter(year >= 1998) %>%
  filter(country_name %in% country_list) %>%
  group_by(aut_ep_id) %>%
  summarize(country_name = first(country_name), 
            start_year = first(year), 
            end_year = last(year), 
            aut_ep_outcome_agg = last(aut_ep_outcome_agg),
            AFI_decline = first(v2xca_academ_lag) - last(v2xca_academ), 
            decline_episode = any(decline_episode ==1)) %>%
  drop_na(aut_ep_id) %>%
  dplyr::select(-aut_ep_id) %>%
  mutate(aut_ep_outcome_agg = case_when(aut_ep_outcome_agg == 1 ~ "Democratic breakdown", 
                                        aut_ep_outcome_agg == 2 ~ "No democratic breakdown", 
                                        aut_ep_outcome_agg == 3 ~ "Regressed autocracy",
                                        aut_ep_outcome_agg == 4 ~ "Outcome censored"))

democracies_autocratiaztion_cases$start_year <- as.character(democracies_autocratiaztion_cases$start_year)
democracies_autocratiaztion_cases$end_year <- as.character(democracies_autocratiaztion_cases$end_year)

democracies_autocratiaztion_cases <- democracies_autocratiaztion_cases %>%
  filter(!aut_ep_outcome_agg == "Regressed autocracy") %>%
  arrange(decline_episode, aut_ep_outcome_agg)

table(democracies_autocratiaztion_cases$aut_ep_outcome_agg)


datasummary_df(democracies_autocratiaztion_cases, output = "outputs/Table_AX2.docx")

## Table 2 ##
  
democracies_autocratiaztion_cases %>%
  group_by(decline_episode, aut_ep_outcome_agg) %>%
  summarize(no = n())





  
  


